home *** CD-ROM | disk | FTP | other *** search
/ MacFormat 1994 August / August CD.bin / Shareware / Education / numericalmethods Folder / chap_9 / hamming.m < prev    next >
Encoding:
Text File  |  1994-06-05  |  899 b   |  34 lines  |  [MATF/MATL]
  1. function [T,Y] = hamming(f,T,Y)
  2. % [T,Y] = hamming(f,T,Y)
  3. % Hamming's solution for y' = f(t,y) with y(a) = ya.
  4. % f  is the function, input.
  5. % T  is the vector of abscissas, input.
  6. % Y  is the vector of ordinates, input.
  7. % Remark. The first four coordinates of T and Y
  8. % must have starting values obtained with  RK4.m
  9. % T  is the vector of abscissas, output.
  10. % Y  is the vector of ordinates, output.
  11. n = length(T);
  12. if n<5, break, end;
  13. f0 = feval(f,T(1),Y(1));
  14. f1 = feval(f,T(2),Y(2));
  15. f2 = feval(f,T(3),Y(3));
  16. f3 = feval(f,T(4),Y(4));
  17. h  = T(2)-T(1);
  18. a  = T(1);
  19. pold = 0;
  20. cold = 0;
  21. for k = 4:n-1,
  22.   pnew = Y(k-3) + 4*h*(2*f1 - f2 + 2*f3)/3;
  23.   pmod = pnew + 112*(cold-pold)/121;
  24.   T(k+1) = a + h*k;
  25.   f4 = feval(f,T(k+1),pmod);
  26.   cnew = (9*Y(k) - Y(k-2) + 3*h*(-f2+2*f3+f4))/8;
  27.   Y(k+1) = cnew + 9*(pnew-cnew)/121;
  28.   pold = pnew;
  29.   cold = cnew;
  30.   f1 = f2;
  31.   f2 = f3;
  32.   f3 = feval(f,T(k+1),Y(k+1));
  33. end
  34.